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Models for the thermal and chemical evolution of the Moon 

We have carried out several studies that explore explanations for the role of chemical 
density variations in Moon’s evolution. Meaningful models for the evolution of the Moon must 
explain a number of important magmatic characteristics. Volcanic activity subsequent to the 
formation of its anorthositic crust was dominated by the eruption of mare basalt. 1) The main 
phase of mare volcanism began -500 Myr after the crystallization of the anorthositic crust and 
continued for -lGyr. 2) The picitic glasses, considered to be representative of mare basalt least 
affected by low pressure, near-surface fractionation, were generated by melting, at 400-600 km 
depth, of a source containing components that, on the basis of the magma ocean hypothesis, 
should have crystallized at much shallower depth during fractionation of the anorthositic crust. 3) 
Mare basalts occur primarily in one region of the Moon. Recent topographic data demonstrate 
that the earlier idea that mare basalt flooded areas of low elevation is not correct. Large areas of 
very low elevation do not contain mare basalt. The hemispheric asymmetry of mare basalt 
distribution on the lunar surface must be explained in some other way (Zhong, et al. 2000). 4) A 
region of the surface roughly correlating with that containing mare basalts also is thought to 
contain high subsurface concentrations of KREEP which was excavated during the formation of 
large impact basins. This so-called Procellarum KREEP Terrane (PKT) is responsible for the 
Imbrium basin-centered thorium anomaly mapped by Lunar Prospector. 

A model in which dense ilmenite and pyroxene-rich cumulates (IC), that crystallize 
during the solidification of the last 5-10% of a magma ocean, sink into the underlying mantle 
may explain important aspects of lunar evolution. The IC may carry with it high 
concentrations of incompatible elements, including heat producing U, Th, and K, that would 
have become progressive enriched in the residual liquid as the magma ocean crystallized. We 
suggested that dense IC cumulates may have differentiated to form a core or a layer 
surrounding an existing metallic core. It is interesting to note that a recent inversion of lunar 
free oscillation data indicates a lunar core with a density appropriate for IC, rather than metal. 
Subsequent thermal expansion due to radiogenic heating would make this material buoyant, 
giving rise to upwelling and melting (Zhong, et al. 2000). This model can thus potentially 
explain the timing of mare volcanism and, if this upwelling developed at spherical harmonic 
degree-one, why mare volcanism in concentrated in one large area of the surface. 


KREEP-rich materials excavated during the formation of the Imbrium basin must have 
pre-existed mare basalts that fill the basin. Both the concentration of KREEP-rich material 
beneath the PKT and later mare volcanism in the same region due to the upwelling and 
melting of initially deep IC cumulates are not a natural consequence of a model of buoyant 
upwelling of a core or deep IC layer unless some preconditioning of the mantle beneath this 
region determined the location of later upwelling. 

It is thus important to explore mechanisms which may regionally concentrate a KREEP 
layer early in lunar evolution and which -400-500 Myr later will give rise to mare basalt 
generation beneath the same region. It has been suggested that melting due to radioactive 
heating of a thickened KREEP layer beneath the PKT may give rise to mare volcanism in that 
region. We explored (Parmentier, et al., 2002) whether gravitational differentiation of a layer 
containing IC cumulates is a possible mechanism to explain the concentration of KREEP 
beneath the PKT. Thickening of the KREEP layer beneath the PKT during differentiation of 
the IC cumulates and the later partial melting of these cumulates due to radiogenic heating 
could have generated the mare basalts that were localized in the region of the PKT. In this 
case, differentated IC cumulates may never have reached depths greater than that of the 
deepest mare basalt source. However, for this explanation to be valid, differentiation of IC 
cumulates in a relatively thin layer beneath the anorthositic crust must have occurred at a 
long wavelength (spherical harmonic degree-one); and the rate of differentiation and the rate 
radiogenic heating must be of an appropriate magnitude to explain the timing and depth of 
origin of mare basalts. We also showed (Hess and Parmentier, 2001) that a KREEP layer 
formed over the whole Moon concentrated at shallow depth beneath the PKT would have 
thermal and magmatic consequences that are not observed. 

To understand the conditions for which this is possible, we examined the Rayleigh- 
Taylor instability of a dense layer near the surface of a sphere, representing the lunar mantle, in 
which the viscosity varies with depth (Parmentier, et al, 2002). Our results show that spherical 
harmonic degree-one is the fastest growing wavelength of instability if the viscosity of the dense 
layer is sufficiently low relative to that of the deeper mantle or if the viscosity of the mantle 
increases with depth. However, the viscosity increase cannot be distributed over a depth that is 
too large. We showed that with a temperature dependent viscosity, that the depth over which 
radioactivity in the IC will heat the mantle, thus reducing the viscosity, is consistent with the 
timing of mare basalt eruption. 

Convective cooling of a compositionally stratified viscous fluid 

Convective cooling of a stably stratified viscous fluid with strongly temperature 
dependent viscosity is important in understanding the thermal and chemical evolution of planets. 
Several different mechanisms may lead to an initially stably stratified mantle. The overturn of 
initially gravitationally unstable lunar magma ocean cumulates could lead to a stably, 
compositionally stratified mantle (Hess and Parmentier, 1995). This could result in a continuous 
Mg/Fe and density stratification over a depth comparable to that of the original magma ocean 
and/or to a more sharply defined layer of very dense, highly radiogenic late stage ilmenite-rich 
cumulates. Alley and Parmentier (1998) describe numerical experiments on the heating from 
below of an initially stable compositional stratification that provides a basis to test the role of 
compositional stratification in the evolution of the Moon. Zhong, et al. (2000) examined a layer 
surrounding a metallic core or a central core of dense, heating producing element-rich silicate 
material as discussed above. 

The formation of crust by decompression melting in upwelling mantle would also result 
in a stably stratified mantle. The amount of melting in an upwelling parcel of mantle increases 
the height over which it decompresses, and increasing depletion of melt decreases the density of 
the residual mantle. Therefore, the progressive accumulation of buoyant mantle, depleted as a 


consequence of crustal genesis by partial melting, should result in stable stratification of mantle 
beneath the lithosphere. We have previous discussed the possible role of such a buoyant depleted 
mantle layer in Venus (Parmentier and Hess, 1992). Recent studies have pointed out that the 
viscosity increase due to depleting small amounts of intergranular water along with melt may be 
larger than originally thought, perhaps as large as a factor of 300. 

For Mars, SNC meterorites record an early magmatic differentiation. Even though Mars 
is a larger planet than the Moon and therefore expected to have a longer, more active convective 
evolution, isotopic heterogeneity recorded in SNCs is as large as those that have been found on 
Moon. Rb-Sr and Sm-Nd isotopes preserve a record of a major depletion event in the mantle, 
indicating mantle has resisted homogenization. Lu-Hf and rare earth element distributions reflect 
small degrees of melting of a source with residual garnet, corresponding to a depth in Mars of at 
least 250 km. If Martian crust, making up approximately 4% of total planetary volume, were all 
generated by adiabatic decompression melting, a significant compositional stratification of the top 
of the mantle would result. This amount of crust would require about 20% melting of the upper 
200 km of Mars. More likely lower degrees of melt distributed over greater depths. How much 
modal cpx and gamet are represented in a chondritic material? 20% melting of garnet or spinel 
bearing mantle would create a density change equivalent to a temperature change of almost 
500°C. Could inhibit convective heat transfer in the mantle. As discussed below, a low mantle 
heat flux may help explain relatively low temperatures at the base of the crust required to explain 
large elastic thicknesses and the long-term preservation of crustal thickness variations 
(Parmentier and Zuber, 2001 ; 2002). It is thus plausible that planets are initially stably stratified 
over substantial depths and that the density variations are large. 

We have carried out numerical experiments examining the instability and subequent 
mixing in a stably stratified fluid layer cooled from above (Zaranek and Parmentier, 2003). The 
obective is to use numerical experiments to derive scaling relationships that can be applied to 
planetary evolution - not detailed models of a particular planet or evolutionary scenario. We 
consider a 2D cartesian domain containing a Boussinesq fluid with infinite Prandtl number. On 
the top, no-slip boundary of this domain, temperature is set at zero. The free-slip bottom 
boundary is adiabatic and the solutions are horizontially periodic. With no volumetric heating or 
heating from below, cooling of the fluid, initially at a uniform temperature 7), creates convective 
instability that is opposed by a stable compositional density stratification with a linear gradient 
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Composition advects but does not diffuse. For thermally activated creep with a large activation 
energy Q, viscosity increases sharply with decreasing temperature thus significantly reducing the 
thermal buoyancy in the boundary layer at the top of the fluid available to drive convective 
motions. We consider an exponentially temperature dependent viscosity 

M(T) = p i exp(e(T i -T)!T t ) 

With e = Q/RT 0 where Q is the activation energy, this relationship provides a good 
approximation to the Arrhenius law with if the temperature range is not too large. Since viscosity 
variations as large as indicated by creep activation energies are difficult or impossible to treat 
with high numerical resolution, the scaling laws that we derive for realizable viscosity variations 
can be extraplolated to stronger temperature dependence. 

Convective instability begins with the initial penetration of cold plumes to a depth 
determined by the compositional stratification. Subsequently formed plumes penetrate to greater 
depth resulting in a mixed layer that cools and thickens with time. The cooling and thickening of 
mixed layer makes it less buoyant than fluid at the top of the underlying, still hot and stratified 


fluid. Continuity of density across the boundary between the bottom of the mixed layer and the 
top the underlying fluid gives 
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where z mixed and z L are the mixed layer and conductive lid thicknesses respectively, T m is the mean 
temperature in the mixed layer, and a is the thermal expansion of the fluid. This relationship 
satisfied well by numerical experiments over a range of conditions, indicating that the mixed 
layer thickens mainly by buoyant, rather than purely viscous, entrainment. 

Cooling of the convecting mixed layer is determined by heat conduction through the lid 
and by the heatflux delivered convectively to its base 
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This relationship follows from earlier work (e.g. Grasset and Parmentier, 1998) where 
AT c =2.23Ti/e\s the temperature over which the viscosity increases by a factor of ten. Cooling of 
the mixed layer, assuming that it thickens fast enough that no heat is conducted into it from the 
still stratified region beneath, is given by the energy balance 
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Combining these equations to eliminate / and T m gives the change in mixed layer thickness with 
time 
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The mixed layer thickens approximately as the squareroot of time so that the expression on the 
right side of the above equation can be considered an effective diffusivity for mixing k m . 
Numerical experiments agree well with this predicted scaling which accounts very well for the 
dependence of mixed layer thickness on yand s. We find values of c in the range of 0.130 to 
0.159 consistent with that for laboratory and numerical experiments with no compositional 
stratification. Some of our numerical experiments also indicate that the mixed layer stops 
thickening, possibly due to the development of a separately convecting layer beneath it. We are 
continuing to study the mechanisms and scaling of the behavior. 

Mars thermal evolution 

Variations in the thickness of a uniform density crust cause lateral pressure gradients that 
generate flow in a ductile lower crust from areas of larger to smaller crustal thickness, thus 
evening out the thickness of the crust over time. Analysis of gravity and topography data for 
Mars provide estimates of crustal thickness showing long wavelength lateral variations present 
even beneath features of great age. The longest wavelength present is the north-south crustal 
thickness variation, but it is not clear that some of the apparent crustal thickness variation 
(derived assuming a constant density crust) does not reflect variations in crustal density. 
Alternatively, the Tharsis Rise appears to be a large center of basaltic volcanism, that began 
forming early in the evolution of Mars, with topography supported largely by thickened crust, 
presumably of relatively uniform composition. Other features such as large basins, like Hellas, 
while apparently younger than Tharsis, preserve shorter wavelength crustal thickness variations. 



In this study we are exploring the implications that the survival of crustal thickness variations 
would have for the thermal evolution of Mars. 

Previous studies have noted that the preservation of crustal thickness variations places 
contraints on either the thickness of the crust through the temperature dependence of crustal 
viscosity and previous studies have investigated constraints on crustal thickness using steady- 
state thermal evolution models which assume a heatflux at the base of the crust equal to the 
radiogenic heat production in the mantle. With a dry diabase flowlaw, crustal thickness as large 
as 100 km will allow the preservation of crustal thickness variations. 

Since the rheology of the lower crust is likely to very strongly temperature dependent, this 
places important constraints on thermal evolution. We have formulated thermal models that 
include both radiogenic heating and secular cooling. These models treat the strong temperature 
dependence of flow by thermally activated creep and the presence of conductive lid above a 
convecting interior using the results of Grasset and Parmentier (1998). Our results (Parmentier 
and Zuber, 2001 ; 2002) indicate that the contribution of the heat flux from secular cooling of the 
core and mantle to temperatures at the base of crust is very significant. For initial temperatures 
for subsolidus cooling in the range of 1500 - 2000°C, and with a crustal thickness as thin as 50 
km and with laboratory flow laws representing a range of plausible crustal compositions, both 
wet and dry, we find that crustal thickness variations even with a lateral scale as large as Tharsis 
cannot be preserved for times on the order of 4 Gyr. We suggest two mechanisms that can reduce 
the temperatures of the lower crust thus preserving crustal thickness variations without requiring 
a very thin crust: compositional stratification of the mantle (Zaranek and Parmentier, 2003) and 
hydrothermal cooling of the crust. 

If the early crust formed by adiabatic decompression melting, for example in mantle plumes, 
the extraction of melt to from crust would leave behind a compositionally buoyant, depleted 
mantle, as discussed above. The results reported above show that stable compositional 
stratification due to crustal genesis may significantly inhibit convective instability at the top of 
the mantle. An estimate of the relaxation rate of crustal thickness variations assuming only the 
conduction of heat from a non-convecting mantle is shown in Figure 4a. Relaxation rates are 
calculated for a dry diabase rheology following the approach outlined by Parmentier and Zuber 
(2001, 2002). Only low values of the fraction of chondritic heat sources fractionated into the 
crust would give relaxation rates smaller than 10' 17 /sec as required to preserve crustal thickness 
variations over the last *4 Gyr. This appears possible only if the crustal heating fraction is less 
that about 20%. Noachian elastic lithosphere thicknesses as large as 50 km (McGovern, et al. 
2001) would also require a relatively low surface heatflux of about 40-50 mW/m 2 . 

As an alternative to low mantle heat flux, a porous, permeable regolith, if water saturated, 
may allow hydrothermal convective cooling of upper crust . On the Moon, low seismic velocity 
in the upper 25 km of the crust is though to reflect fractured, porous rock with a porosity of about 
1%. An impact-fractured and brecciated layer should also have been present on Mars perhaps 
with a similar porosity to that on the Moon. The bottom of such a porous layer is determined by 
the close of fractures under hydrostatic stress, then the layer thickness on Mars should be about 
1/3 that on the Moon. Hydrothermal cooling is determined a Rayleigh number that depends on 
the product of the permeable and layer thickness and the physical properties of water. Assuming 
that fractures occur a planar orthogonal network of cracks, the permeability is determined by the 
product of the porosity and the square of the crack spacing. For fractures spaced 1 cm apart, 
perhaps typical of the clast sizes in lunar breccia samples, and using the appropriate physical 
properties of water, values of Nu (Nusselt number), the ratio hydrothermal convective heat 
transfer to purely conductive heat flux , would about a factor of ten. For Nu=10 and an 8 km 
layer thickness, the relaxation time for a mantle heatflux expected early in Martian history (60-80 
mW/m 2 ) is small enough to preserve crustal thickness variations. 

Hydrothermal alteration may also explain the difference in the intensity of magnetic 
anomalies between the southern highlands and northern plains. If crustal mineralogy is 



responsible for remnant magnetization and if hydrothermal alteration extends to the same depth in 
both the highlands and the plains, then thicker crust in the highlands would result in a thicker, 
unaltered, still magnetized layer. Hydrothermal ly alterated minerals may also fill cracks reducing 
the permeability. Hydrothermal cooling may thus be most effective when impacts are frequent 
enough to continually generate new fractures. 
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